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ABSTRACT 

We use 3D SPH calculations with higher resolution, as well as with more realistic 
viscosity and sound-speed prescriptions than previous work to examine the eccentric 
instability which underlies the superhump phenomenon in semi-detached binaries. We 
illustrate the importance of the two-armed spiral mode in the generation of super- 
humps. Differential motions in the fluid disc cause converging flows which lead to 
strong spiral shocks once each superhump cycle. The dissipation associated with these 
shocks powers the superhump. We compare 2D and 3D results, and conclude that 3D 
simulations are necessary to faithfully simulate the disc dynamics. We ran our sim- 
ulations for unprecedented durations, so that an eccentric equilibrium is established 
except at high mass ratios where the growth rate of the instability is very low. 

We collate the observed data on superhumps. Our improved simulations give a 
closer match to the observed relationship between superhump period excess and bi- 
nary mass ratio than previous numerical work. The observed black hole X-ray transient 
superhumpers appear to have systematically lower disc precession rates than the cat- 
aclysmic variables. This could be due to higher disc temperatures and thicknesses. No 
high-resolution 3D disc with mass ratio q > 0.24 developed superhumps, in agreement 
with analytical expectations. 

The modulation in total viscous dissipation on the superhump period is over- 
whelmingly from the region of the disc within the 3 : f resonance radius. The pre- 
cession rates of our high resolution 3D discs match the single particle dynamical pre- 
cession rate at 0.87i?3 : i. As the eccentric instability develops, the viscous torques are 
enhanced, and the disc consequently adjusts to a new equilibrium state, as suggested 
in the thermal-tidal instability model. We quantify this enhancement in the viscosity, 
which is ~ 10 per cent for q = 0.08. The disc motions can be described as superpo- 
sitions of the S(k,l) modes, and the disc executes complex standing wave dynamics 
which repeat in the inertial frame on the disc precession period. We characterise the 
eccentricity distributions in our accretion discs, and show that the entire body of the 
disc partakes in the eccentricity. 

Key words: accretion, accretion discs — hydrodynamics — instabilities — methods: 
numerical — binaries: close — novae, cataclysmic variables 



1 INTRODUCTION 

Cataclysmic variables (CVs) are semi-detached binaries with 
a Roche- lobe filling low-mass donor star, mass M%, transfer- 
ring matter onto a white dwarf (WD) primary, mass Mi, 
via an accretion disc. The SU UMa-type dwarf novae (DNe) 
are short-period cataclysmic variables which display two 
distinct modes of outburst. The normal outbursts are at- 
tributed to a thermal-viscous limit-cycle between low and 
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high viscosity states (Osaki 1974; Hoshi 1979; Meyer & 
Meyer-Hofmeister 1981), and consequently between low and 
high mass transfer states (see Lasota 2001, for a review). 
The larger amplitude and longer-lasting superoutbursts are 
characterised by a periodic photometric modulation known 
as superhumps. Superhumps are attributed to an eccentric 
apsidally-precessing accretion disc (Vogt 1982). The super- 
hump period, P s h, is a few percent longer than the orbital 
period, P or h'- the orientation of the mass donor star relative 
to the progradely precessing eccentric disc repeats on P s h. 
Whitehurst (1988) and Lubow (1991a) explained that a disc 
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which encounters a 3:1 eccentric inner Lindblad resonance 
with the tidal potential of the secondary star may become 
eccentric and precess. A mass ratio g = M2/M1 < 1/4 
is required for the tidal truncation radius, i?tidcs, to lie 
outside the 3:1 resonance radius, R3.1 (Paczynski 1977). 
The thermal-tidal instability (TTI) model (Osaki 1989) at- 
tributes the increased brightness and duration of superout- 
bursts over normal outbursts to an enhanced viscous torque 
acting once the disc becomes eccentric. 

Superhumps arise in several guises, summarised by Pat- 
terson et al. (2002a). 'Common' or 'normal' superhumps en- 
sue after the onset of superoutburst in SUUMa systems. 
They are powered by the periodically varying tidal interac- 
tion which modulates dissipation in the disc (e.g. Foulkes 
et al. 2004). 'Late' superhumps sometimes follow common 
superhumps and are roughly anti-phased to them. The mod- 
ulation in energy dissipation at accretion stream's impact on 
the non-axisymmetric disc powers late superhumps (Rolfe, 
Haswell & Patterson 2001). 'Negative' superhumps, with P s h 
slightly less than P or b, are sometimes observed simultane- 
ously with the more usual 'positive' superhumps, and may 
be related to retrograde precession of a warped accretion disc 
(Patterson et al. 1993b; Murray & Armitage 1998; Foulkes, 
Haswell & Murray 2006). 'Permanent' superhumps are seen 
in high mass transfer rate systems: the nova-like variables, 
old novae and some AM CVn systems. Some low mass X-ray 
binaries (LMXBs, the neutron star and black hole analogues 
of CVs) also show superhumps (Bailyn 1992; Haswell 1996; 
O'Donoghue & Charles 1996). In LMXBs optical emission 
arises overwhelmingly from the reprocessing of X-rays, and 
superhumps arise from a modulation in reprocessing caused 
by the changing solid angle subtended by the tidally flexing 
eccentric disc (Haswell et al. 2001). Recently superhumps 
were reported in the microquasar GRS 1915+105, which has 
an orbital period exceeding 30 days (Neil, Bailyn & Cobb 
2006). 

In all the above cases, the fractional superhump excess, 
e = (P s h — Port>)/Porb varies with q, with |e| increasing with 
higher values of q. The exact relationship has proved difficult 
to determine. We have performed unprecedentedly compre- 
hensive numerical simulations of apsidally precessing accre- 
tion discs, and we present them in the context of previous 
numerical work, the observational data and salient analyt- 
ical theory. Section 2 describes our simulations and exam- 
ines the growth rates of the eccentric instability; the en- 
hanced viscous torques; the superhump light curves which 
result from the eccentric instability; and compares 2D and 
3D simulation results. In section 3 we use two methods to 
quantify the eccentricity distributions in our simulated discs. 
Section 4 focuses on the e — q relationship. In section 5 we 
discuss our findings. Section 6 gives a summary list of our 
principle conclusions. 



2 SPH SIMULATIONS OF A PRECESSING 
ACCRETION DISC 

SPH is a Lagrangian method which models fluid flow as a set 
of moving particles. A detailed review is given by Monaghan 
(1992). SPH simulations by Murray (1998) provide consid- 
erable support for the TTI model, showing that the energy 
released from a disc that has become tidally unstable is suffi- 



cient to account for the excess luminosity of a superoutburst. 
Foulkes et al. (2004) carried out 2D SPH simulations of a 
binary system with mass ratio 0.1. They show an eccentric, 
non-axisymmetric precessing disc of changing density, which 
is continuously flexing and relaxing on the superhump pe- 
riod. Very clear too in the surface density maps are tightly 
wrapped spiral density waves which extend from the outer- 
most regions to small radii. They produce shear and dissi- 
pation in the outer disc, and propagate angular momentum 
outwards, allowing disc gas to move inward. 

2.1 Simulation details 

For the calculations presented here, an SPH code is used 
which has been designed specifically for accretion disc prob- 
lems. Detailed description of the code can be found in Mur- 
ray (1996, 1998). The code is normalised in units of a, the 
binary separation, for distance; Aft = M1 + M2 for mass; and 
Porb/(27r) for time. This code has since been updated to in- 
clude adaptive spatial resolution, allowing the SPH smooth- 
ing length A to vary in both space and time (Murray, de Kool 
& Li 1999). Here, A is set to a maximum value of 0.005 a. The 
code has also been extended to three dimensions (Murray & 
Armitage 1998). 

Simulations were run for a range of mass ratios, de- 
tailed in Table [1] The simulations were built up from zero 
mass and a single particle was injected at the L\ point each 
timestep At, so simulating the mass-transfer stream from 
the secondary. We ran simulations at different mass resolu- 
tions, where At is between 0.01 Q~ T \ and 0.0025 H^. This 
latter value is at resolution higher than previous calculations 
(Murray 1996, 1998; Foulkes et al. 2004), and we find our 
results stable to a further reduction in At. The number of 
particles in the accretion disc in these high-resolution simu- 
lations is of order 100 000, and the average number of 'neigh- 
bours', i.e. the average number of particles found within 2A 
of each other that are used in the SPH update equations, is 
between 20 and 30. 

The stream boundary conditions at L\ are a function 
of q, calculated by Lubow & Shu (1975) from perturbation 
analysis of L\. We follow their calculations to determine 
the direction prograde of the binary axis at which particles 
are to be injected, #i n j. The initial speed of the particles, 
^inj = 0.1aSl or b, is determined from c s at Li, where the 
z- velocity is an arbitrary fraction of this. The actual value 
is not critical as the gas will rapidly accelerate to become 
supersonic, and the stream is not well-resolved in our sim- 
ulations. For some simulations the third dimension is sup- 
pressed. For the inner boundary condition a hole of radius 
Ri centred on the primary was used, with Ri set to values 
between 0.03 a and 0.1a, and particles entering this were 
removed from the simulation. 

The calculations presented in Murray (1998) and Mur- 
ray, Warner & Wickramasinghe (2000) were of cool isother- 
mal discs. Here we model a more realistic steady-state disc 
where c s is a function of disc radius, r, and is given by 
c s = Cr" 3 ' 8 . C is a constant and in general we have set 
it equal to O.05af2 or b. This means that c s at the resonance 
radius in each of our simulated discs will be ~ 0.067 a fl or b as 
detailed in Table H The SU UMa systems Z Cha and OY Car 
in outburst have brightness temperatures, Tbr, in the outer 
disc of ~ 6000 - 7000 K (Home & Cook 1985; Bruch, Beele 
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Table 1. Summary of simulations and results. The first column denotes the run number. Columns 2 to 6 describe the simulation 
parameters, recording respectively whether the simulation is conducted in 2D or 3D, the mass ratio, the injection time step, the constant 
describing the sound speed (c s = Cr -3 / 8 ) and the radius of the primary (the central hole). The time at which the simulation was 
terminated is given in column 7. The remaining columns record outcomes of the simulations: the mean superhump excess as measured 
from the simulated lightcurve, the final total number of particles in the simulation, the final average number of neighbours, the strengths 
of the eccentric and 2-armed spiral modes averaged over the final 10 superhump periods in each simulation, and, as a measure of when 
the disc initially encounters the resonance, the time at which S(i,o) = 0.01. f = 1.0 in all cases. 
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Baptista 1996). For the eight SUUMa eclipsing systems 



in Table [3] we obtain 0.031 a 0, o 



< O.036af2 or b at 



the disc midplane at R3.1, assuming a mass transfer rate 
of M — 10 -9 Mq yr _1 and a fully ionised cosmic mixture 
of gases. Hence our simulations have a sound speed within 
a factor of two of that prevailing in reality, a much better 
match than previously possible. 

We have also improved on the viscosity used in previous 
calculations. Our code includes an artificial viscosity term 
which generates a shear viscosity in the disc, 



u(r) = k(c s H, 



(1) 



where £ is the dimensionless artificial viscosity parameter. 
C = 1 here, and H is the disc scaleheight. k may be found 
analytically, and for a standard cubic spline kernel in three 
dimensions, « = 1/10, whilst in two dimensions, k = 1/8. 
The bulk viscosity is fixed to be twice the shear viscosity. 
Using the Shakura-Sunyaev viscosity parametrisation (y = 
aCsH) (Shakura & Sunyaev 1973), our simulated 3D and 
2D discs have a(i?3,i) = 0.1 and 0.125 respectively. This 



matches estimates of a ~ 0.1 — 0.2 derived from observation 
of systems in the high viscosity state (Smak 1999). 



2.2 Simulation Results 

Table Q] summarises our simulations. In some cases these 
have run for > 2000 orbits without reaching mass-transfer 
equilibrium; one has run for over a year of elapsed time. 
The rate of energy dissipation (i.e. viscously-generated lu- 
minosity with units Ma 2 fio rb ) from different regions in the 
disc is recorded each timestep and used to produce simula- 
tion lightcurves. P s h was determined from timings of maxima 
in dissipation in a smoothed lightcurve for the disc region 
r > 0.3 a. Column 8 of Table \T\ gives the mean value of e 
so obtained, over the time in which the system has reached 
equilibrium where applicable. Figure Q] shows the disc evo- 
lution for the 3D simulations 1 to 12. Here the disc mass is 
taken to be the number of SPH particles in the disc. The disc 
eccentricity is estimated from the eccentric mode strength, 
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Table 2. Characterisation of the simulated accretions discs at the 3:1 resonance radius. The 3:1 resonance radius is recorded in the 
fourth column, followed by the sound speed, the scale height, the characteristic value of the ratio of the disc semi-thickness to the radius 
as given in Goodchild & Ogilvie (2006), the shear viscosity, the bulk viscosity, the Shakura-Sunyaev parameter for the shear viscosity 
and that for the bulk viscosity. 
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2.2.1 The growth of the eccentricity 

Figure [T] shows that initially, the eccentricity grows expo- 
nentially. Lubow (1991a) found an exponential eccentricity 
growth rate which is proportional to the square of the mass 
ratio; in contrast we see a very low growth rate for high 
mass ratios (g > 0.2) (Figure [2]). Our results can be prof- 
itably compared with those of Goodchild & Ogilvie (2006) 
who formulated a single equation to describe the resonant 
excitation, propagation and viscous damping of the eccen- 
tricity in a 2D accretion disc. They showed that the res- 
onance may have the effect of locally suppressing the ec- 
centricity which consequently leads to extremely low eccen- 
tricity growth rates. We find for q = 0.2195, for example, 
a steady state is not reached until ~ 1700 orbital periods, 
and for q = 0.2346 a steady state is still not achieved after 
~ 2500 orbital periods. 
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Figure 2. Eccentricity growth rate (growth rate of the strength 
of the (1,0) eccentric mode) as a function of mass ratio for sim- 
ulations presented in this work, and in previous works (Murray 
1998). Simulation parameters are as indicated in the legend. 
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Table 3. Parameters of the eight eclipsing SUUMa systems: mass ratio, primary mass, total system mass (Mi + M2), observed disc 
precession rate (from the measured e value), dynamical precession rate as calculated from Equation[6]at the location of the 3:1 resonance, 
inferred pressure contribution to precession (w bs — w dyn) an d pitch angle of the spiral wave (see Section 4 for details). Errors on u> oos 
and i p correspond to the range of superhump periods observed. In the last but one column the sound speed, in SPH units, is calculated 
for the midplane of each disc at the 3:1 resonance radius, assuming a mass transfer rate of 10 — 9 Mq yr~ 1 and a fully ionised 'cosmic' 
mixture. References for Mi and Mt are provided in the final column. Details and references for all other observational data can be found 
in Table [5] 



System 


q 




Mi 


M t 


^oba 

(radd^ 1 ) 


Wdyn(^3:l) 

(radd- 1 ) 


LJpv 

(radd- 1 ) 


ip 
(°) 


C B (R 3 :l) 

afl orb 


Ref 


OY Car 





102 ± 0.003 


0.685 ±0.011 


0.755 ±0.011 


2.189t° £I 


3.649 


-1.460 


12 69+ 2 ' 44 


0.0334 
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XZ Eri 





1098 ± 0.0017 


0.767 ±0.018 


0.851 ±0.018 


AOSI -0.0035 


4.015 


-1.319 


ro.uz_ 16 


0.0321 
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IY UMa 





125 ± 0.008 


0.79 ± 0.04 


0.89 ±0.04 
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-1.506 


11 17+ 1 - 42 


0.0322 
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1495 ± 0.0035 
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3.161±° S34 


4.282 
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HT Cas 





15 ±0.03 


0.61 ± 0.04 
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4.343 


-1.618 


11.62^° 


0.0347 
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9 34Q+0-080 
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-1.389 
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±U.D<i_ g 2 
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2 


OU Vir 





175 ±0.025 


0.90 ±0.19 


1.06 ±0.19 


2.732+0.033 


4.987 


-2.255 


8.76t°' 06 


0.0305 
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V2051 Oph 





19 ±0.03 


0.78 ± 0.06 


0.93 ±0.07 




6.201 


-3.453 


7 oo + l.H 
' • oo -0.33 


0.0312 


7 



1 Wood et al. (1989) 2 Feline et al. (2004c) 3 Steeghs et al. (2003) 4 Wade & Home (1988) 5 Home, Wood & Stiening (1991) 
6 Feline et al. (2004a,b) 7 Baptista et al. (1998) 



2.2.2 The enhanced viscous torques of an eccentric disc 

For all calculations in which the disc (eventually) becomes 
sufficiently eccentric we see similar behaviour but on dif- 
ferent timescales. As Figure [T] shows, initially the disc mass 
builds, exponentially approaching a steady-state value. Then 
as eccentricity increases, the disc mass tends to a new lower 
steady state value. This reveals the non-eccentric accretion 
disc approaching an equilibrium between the tidal removal 
of angular momentum and the angular momentum added 
by material from the LI point. As the disc becomes eccen- 
tric it readjusts to a new equilibrium in which tidal removal 
of angular momentum is more efficient. This is exactly the 
premise of the TTI model. 

Table U quantitatively examines this. As a measure of 
the increase in efficiency of tidal removal of angular momen- 
tum in the eccentric disc, we took the decrease in disc mass 
between maximum and the value it finally reached in equi- 
librium (column 4 of Table [4}. For a steady state accretion 
disc the total mass is 



Mtot = 



/' 

JR. 



2RM 
MB) 



R, 

Rout 



dR 



(2) 



(Frank, King & Raine 1985). The mass transfer rate through 
Li remains constant, and we are comparing the equilibria 
with and without an eccentric precessing disc. In both cases, 
the mass accretion at all disc radii must equal the mass 
transfer rate at Li. Thus, the change in disc mass can then 
be related to a change in viscosity by 



Mb - A/ a = AAf = A 



/ (R) dR f f (R) dR 



(3) 



MR) J MR) 

where A is a constant, and the subscripts b and a respec- 



tively refer to before and after the disc became fully eccen- 
tric. Defining 



1 



"b(-R) 



ff(R)dR> 



(4) 



where is some unknown weighted mean value of v D (R), 
and similarly for -U, then we have 



AM 



Aff(R)dR - - 

This quantity is given in column 5 of Table [4] The radial de- 
pendence of H is c s (r//i) 1/ ^ 2 r, where jj, = l/(g+l). Together 
with the sound speed prescription that we use, then this 
allows us, via Equation [1] to give explicitly the fractional 
change in viscous torque necessary to bring about the de- 
crease in disc mass we see (column 6 of Table We assume 
Rout — -Rtidos- Here, we use the formulation i?tidcs — 0.9iii 
(Frank, King & Raine 1985) , where Ri is the effective Roche 
lobe radius (Eggleton 1983), though we note the limitation 
of this approximation (e.g. Murray, Warner & Wickramas- 
inghe 2000; Truss 2007). For runs 6, 7, 9 and 10 we measured 
the average outer disc radius, averaging the position of the 
outermost particle as a function of azimuth over a super- 
hump period. We found in each case the difference between 
this value and -Rtides to be < 0.006 a. 

2.2.3 The superhump 

The development of the superhump for two selected simula- 
tions is shown in Figure[3] For q — 0.0526, superhumps were 
only temporarily present at a time when the eccentricity was 
highest. For q = 0.2121 we see that the superhump profile 
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Figure 1. Disc mass (red) and disc eccentricity (green) as a function of time for the 3D simulations 1 to 12. Every 5000th and 2000th 
timestep is plotted respectively. Also plotted is the superhump period (blue) in the form of O-C. This is calculated relative to the mean 
super hump period from each run, derived from the point at which the superhump signal becomes well-formed. The dissipation (smoothed) 
from the disc for r > 0.3 a is shown in grey. The mass ratio is shown in the upper left of each panel. 



evolves with time as the disc is reaching eccentric equilib- 
rium. In the case of other more extreme mass ratios (smaller 
values of q) there are higher frequency periodic components 
present in the early development of the superhumps (Smith, 
PhD thesis in prep.). Our discs accumulate from zero mass 
with no switch between viscosity states, so the development 
of their superhump will differ from that of discs observed in 
superout burst. 

Variations in the superhump period are displayed in the 
form of O-C ('observed' minus 'calculated') in Figure [T] In- 
tervals of non-zero 2nd derivative in the value of O-C indi- 
cate period changes. Generally P s h decreases as the system 
approaches eccentric steady-state. For q — 0.0526 a steady- 
state eccentric disc is not achieved, and the period evolution 
differs. 

Each simulated lightcurve was folded on the mean su- 
perhump period once equilibrium was reached (Figure [4}. 
These are asymmetric with, for most, a steep rise and slower 
decline. A secondary hump structure is seen, the profile dif- 
ferent for different mass ratios. In general, these superhump 
profiles resemble those observed in CVs cf. figure 8 of Patter- 
son et al. (1995), figure 4 of Imada et al. (2006) and figure 6 
of Maehara, Hachisu & Nakajima (2006). 

In Figure [5] we show the density distributions of the 



accretion disc for run 9 at selected superhump phases, and 
compare these with the density profile of a disc which does 
not show superhumps. The disc is not maximally distorted 
at the time of maximum energy production (0 s h = 0) as 
might be expected, but at <f> s h = 0.71. For q = 0.1111, this 
coincides with a small secondary peak seen in Figure [4] The 
superhumping disc is most similar to the non-superhumping 
disc at superhump maximum, the most visible difference be- 
ing in the strongly enhanced spiral density waves in the su- 
perhumping disc. Furthermore, the appearance of the spiral 
density wave changes dramatically over P s h- This illustrates 
the crucial importance of the spiral density waves to the 
superhump phenomenon, as suggested by Lubow (1991a) 
and Osaki (2003) who considered analytic theory and obser- 
vations respectively. Comparing superhump maximum with 
superhump minimum, then we see more a more open spiral 
structure at superhump maximum with outer reaches show- 
ing enhanced density. After superhump maximum, differen- 
tial motion in the superhumping disc causes the spiral to 
become eccentric and more tightly wrapped with lower den- 
sity contrast (see the bottom two panels in Figure [5}. Each 
radius in the disc has its own characteristic Keplerian angu- 
lar velocity, and eccentric mode precession rate. Differential 
precession in the eccentric fluid disc cannot persist, because 
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Table 4. Enhanced viscous torque resulting from the accretion 
disc encounter with the eccentric resonance. The fourth column 
records the fractional decrease in disc mass before and after the 
disc becomes fully eccentric. The fifth column gives a measure of 
the change in viscous torque indicated by this mass decrease as 
described in the text and given by Equation [5] The final column 
assumes the radial dependence of viscosity is given by Equation[T] 
and records, in that case, the fractional change in viscous torque 
necessary to bring about the decrease in disc mass. 
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Figure 4. Equilibrium lightcurves (energy dissipation rate for 
radii r > 0.3 a) folded on the derived mean superhump period for 
runs 5 to 10. The superhump cycle is repeated for clarity. 
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Figure 3. Evolution of the simulated lightcurve for runs 6 and 
11. Each panel covers 5 orbits. The lightcurve, which has been 
smoothed somewhat, is given by the rate of energy dissipation 
for radii r > 0.3 a. The mass ratio is shown in the uppermost 
panel of each. 



this would cause widespread orbit crossings. This converg- 
ing fluid motion causes the strong spiral shock shown in the 
upper right panel of Figure[5] and it is the dissipation associ- 
ated with this shock which powers the observed superhump. 
Figure [S] illustrates the mechanism by which spiral waves 
power superhumps. 




Figure 5. Density profiles comparing the non-superhumping disc 
of run 2 (upper left panel) with the supcrhumping disc of run 9 
at 3 selected superhump phases, ^j,. The coordinates are centred 
on the binary system centre of mass. The solid line is the primary 
Roche lobe and the dashed line is the 3:1 resonance radius. The 
upper left panel also shows our definition of azimuth used in Sec- 
tion [37T] The secondary star is at an azimuth of n radians with 
respect to the primary, the position of which is marked with a 
cross. 
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Figure 6. Comparison of results for 2D and 3D simulations with q = 0.1111, runs 24 and 9 respectively. The following are shown above 
and below respectively for each: the disc evolution, details of the simulated lightcurves, folded superhump lightcurves and density profiles. 
Here symbols are as in Figure [5] and the colour scale used is the same for both the 2D and 3D discs. The simulated lightcurves cover 5 
orbital cycles once the disc mass and eccentricity had stabilised. The raw lightcurve is shown in the top panel and increasingly smoothed 
lightcurves below. 



2.2.4 2D versus 3D 

In Figure [6] we compare the results of q = 0.1111 calcu- 
lations in 2D and 3D. In the 2D simulations the accretion 
disc achieves a far higher eccentricity more quickly. This 
in turn affects the superhump profile. We see from Table [1] 
that these 2D accretion discs (runs 19 to 27) can become 
eccentric at much larger values of mass ratio than in the 3D 
simulations, and well beyond the bounds that theory sug- 
gests, indicating that three dimensions and the high mass 
resolution (At — 0.0025) is necessary to accurately simu- 
late the processes involved. In this case, (3D, At — 0.0025) 
no disc develops superhumps where q > 0.24. We also see 
that for the same value of q, the 2D discs precess at smaller 
rates than the 3D discs; compare for example run 24 with 
run 7. 



3 ECCENTRICITY DISTRIBUTION IN 
ACCRETION DISCS 

An eccentric accretion disc underlies the superhump phe- 
nomenon, but how eccentric does the disc become? How does 
the eccentricity vary throughout the disc, and how does it 
change with superhump phase? 

We plotted an estimate of the disc eccentricity in our 
simulations in Figure [T] Calculating the strength of the 
(1,0) mode takes into consideration the disc as a whole. 
For those 3D discs which reached equilibrium and showed 
superhumps, final values of eccentricity e between ~ 0.09 
and ~ 0.23 were seen, with more extreme mass ratios har- 
bouring more eccentric discs. This agrees with modelling 
of line profiles in AM CVn assuming a constant eccentric- 
ity throughout the disc: Patterson, Halpern & Shambrook 
(1993c) found e = 0.1 — 0.2. In their analytic study Good- 
child & Ogilvie (2006) examined the spatial distribution of 
eccentricity, but explicitly avoided examination of the be- 
haviour on the orbital timescale. Their eccentricity distri- 



bution was locally suppressed by the presence of dynamical 
resonances. Next we examine the eccentricity distributions 
of our simulated discs, using two complementary methods 
to characterise their spatial and temporal variation. 

3.1 Eccentricity distribution from particle 
trajectories 

An instantaneous snapshot of the radial eccentricity distri- 
bution was found by projecting the elliptical orbit of each 
particle in the disc using its position and velocity. Particles 
in the mass transfer stream were discounted. Each parti- 
cle was assigned a radius given by an average along this 
elliptical path, weighted by the time spent at each radius. 
This method calculates the trajectories particles would have 
if they orbited an isolated primary star. This simplification 
may introduce certain artefacts into the calculated eccentric- 
ity distributions, but does give an adequate approximation. 

Figure [7] shows the results for each of the 3D simula- 
tions 1 to 12, at superhump maximum where applicable. In 
this and subsequent similar figures, each point represents 
the mean eccentricity of particles, obtained in the manner 
described above, in each of a number of azimuthal and ra- 
dial bins of size 0.l7rrad and 0.1a respectively. Overplot- 
ted on these figures are the 3:1, 4:1 and 5:1 resonance radii 
from right to left respectively. For simulations 11 and 12 the 
2:1 resonance is also shown. Each particle is colour coded 
according to the azimuthal position of the particles it rep- 
resents, as shown in the key, where 9 is the angle defined 
in Figure [S] Purple and blue particles are approaching the 
mass donor star, while green and yellow particles are re- 
ceding from it. The eccentricity distributions at different 
azimuths are clearly distinct. Caution is required, however, 
in interpreting the radii assigned in this section. The parti- 
cles at r > i?3 : i are concentrated at azimuths close to the 
x-axis, where in fact the disc edge is relatively close to the 
compact object: these particles have larger radii assigned to 
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Figure 7. Instantaneous eccentricity distribution plotted as a function of radius for 3D runs 1 to 12, once equilibrium is reached (see 
text for details). Each point represents the average of all disc particles in each of 20 equal azimuthal sections and binned in radius. These 
are colour coded by their azimuthal position in the disc, as defined in Figure [5] according to the scale on the right. The time is shown 
in the upper right of each panel, together with the superhump phase where applicable. 



them than those found further out at other azimuths be- 
cause they have relatively circular projected orbits. 

For all the discs with q ^ 0.0811 the eccentricity in the 
outer disc has a wide spread around e ~ 0.2. This agrees 
reasonably with the eccentricities found in the outer discs of 
OYCar: Hessman et al. (1992) found e = 0.38±0.10; and IY 
UMa: Patterson et al. (2000a) estimate e = 0.29 ± 0.06. For 

0. 2121 ^ q ^ 0.3333 the eccentricity distributions are low 
for the disc within the 5:1 resonance radius, then increase 
at larger radii to e up to ~ 0.3, this apparent eccentric- 
ity likely due to tidal distortions in the outer disc. For this 
range of q, if every particle were individually plotted, the ec- 
centricity distribution outside R5.1 is "double horned", with 
a concentration of points near the upper and lower edges 
of the distribution and a dearth of points midway between 
the envelope. These two "horns" correspond to the halves of 
the disc which are approaching and receding from the mass 
donor star. 

Figure [7] demonstrates that the body of the disc within 
the 5:1 resonance appears eccentric for 0.081 ^ q ^ 0.2195. 
For these discs the entire mass appears to partake in the ec- 
centric instability. This roughly flat region of the eccentricity 
distribution has an eccentricity dependent on q, increasing 
to ~ 0.12 as we move to smaller q. 

Figure [H] shows how the eccentricity distribution devel- 
ops in run 9 as the disc becomes eccentric. Before super- 
humps set in (top left panel) the eccentricity distribution 
resembles those seen in the less extreme mass ratio systems 

1. e. the top row of panels in Figure [7] The remaining 5 pan- 
els of Figure [5] span the interval where the disc eccentricity 
is growing and the disc is emptying as it approaches the 
new mass transfer equilibrium caused by the enhanced vis- 




rcidius 



Figure 8. Development of accretion disc radial eccentricity dis- 
tribution for run 9. Azimuthal colour-coding is as in Figure 171 

cous torque (c.f. Figure [T}. In each case superhump max- 
imum is plotted. Within _Rs : i the eccentricity of the body 
of the disc grows, while the range of eccentricities at any 
given radius remains more or less unchanged. Outside R5.1 
the eccentricity distribution becomes more scattered, as the 
range of eccentricities at any given radius increases. The last 
panel of Figure [S] shows the distribution when the strength 
of the (1,0) mode approaches maximum, but mass equilib- 
rium is not yet reached. It resembles the q = 0.1111 panel 
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Figure 9. As Figure [7] showing the changing radial eccentricity 
distribution with superhump phase, for run 9. 



in Figure [7] which occurs almost 230 orbits later, after mass 
equilibrium is established. 

The eccentricity distribution is shown as a function of 
superhump phase in Figure [9] At some phases and radii 
there is little spread in e, while at other phases there is a 
large spread in e at the same radius. This behaviour may 
be related to the cusps at the resonance points discussed by 
Goodchild & Ogilvie (2006), but their treatment explicitly 
excluded the behaviour of the disc on timescales comparable 
to or shorter than the orbital period. Numerical simulation 
facilitates examination of the flexing of the disc during a 
single orbit, which is important as it is this flexing which 
gives rise to the modulation in viscously-generated or repro- 
cessed light which constitutes an observed superhump. The 
flexing of the disc over the superhump period appears most 
dramatic at radii between and R3.1, in accordance with 
Pearson's (2006) point that the dynamic precession rate at 
the 4:1 resonance actually agrees far better with the obser- 
vations. 

Figure [5] neatly illustrates the changes which occur as 
the binary's gravitational potential moves relative to the 
disc. Even as far in as r=0.15a we see the more eccentric 
particles belonging to the approaching side of the disc at 
superhump maximum and the receding part of the disc at 
superhump minimum. 



3.2 Eccentricity from the mass distribution 

How good a representation of the disc eccentricity are the 
distributions we calculated in section 13.11 ? If we were to fol- 
low a single particle as it orbits in the accretion disc it would 
not move on the elliptical orbit we extrapolated from its in- 
stantaneous velocity and position. We now use an alterna- 
tive way of examining the eccentricity of the accretion disc, 
looking at the mass distribution. 

The disc in run 9 was split into 100 azimuthal sections. 
For each, we recorded the number of particles contained 
within each radial step outwards. In Figure 1101 we show 
the contour maps which result. Every step of 24 particles 
was recorded and colour-coded, so each colour represents a 
'contour' of enclosed mass. To each contour we fitted an el- 
lipse which has one focus at the WD. The fitted parameters 
are the semi-major axis, a som i, the eccentricity, e, and the 





-0.6 -OA -0.2 0.0 0.2 0.4 



-0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 



Figure 11. As Figure [Tol for the non-super humping disc of run 
1. 



angle the semi-major axis makes with the positive s-axis in 
an anti-clockwise direction, ao, measured in radians. This is 
analogous to our definition of 8 in Section [3.11 For the out- 
ermost complete contour, that is the outermost contour for 
which each azimuthal section is represented, marked in dark 
blue, these fitted parameters are displayed in the bottom 
right of each panel in Figure \W\ and apply to the overplot- 
ted magenta ellipse. 

Figure [lOl allows us to look at the disc mass distribution 
in a more quantitative way than simply looking at the den- 
sity distribution. The eye is drawn to the pronounced non- 
axisymmetry outside the 3:1 resonance, though this mass 
constitutes less than 9 per cent of the mass in the disc, and 
is the lowest density region. This mass contributes only 4 
per cent of the total dissipation, and this contribution to 
the modulation on P s h is not in phase with the overall su- 
perhump. The region of the disc within R3.1 is overwhelm- 
ingly responsible for generating the dissipation-powered su- 
perhump in the simulations. 

The four panels in Figure [10] are equally spaced in su- 
perhump phase. It is noticeable that the outermost disc does 
not simply precess as the binary frame moves. Instead the 
flexure of the disc combines with the orbital motion of the 
binary frame to leave the outer edge of the disc almost fixed 
in the binary frame between 4> s h = 0.02 and <f) s h = 0.27; 
similarly the outer edge of the disc remains almost the same 
in the two lower panels at <j) a h = 0.52 and 4>sh = 0.77. Fig- 
ure [11] shows that a non-superhumping disc is extended at 
similar azimuths. These extended disc edges are analogous 
to the raised tides in Earth's oceans. As Figure [TT] shows, 
this effect produces an elliptical shape centred on the pri- 
mary. This illustrates the distinction that needs to be made 
between different contributions to the eccentricity distribu- 
tions found here and in Section [3.11 In the outer disc tidal 
distortions are important, and are present in discs at all 
mass ratios. This is to be distinguished from the (1,0) ec- 
centric mode eccentricity found in the superhumping discs 
which has the primary at one focus. Truss (2007) finds that 
for discs which have not yet become tidally unstable, the 
exact azimuth of the extended 'wing' depends on the mass 
ratio, viscosity parameter and sound speed of the gas (the 
major axis of the outer streamline moves clockwise with de- 
creasing q, or increasing a and c a ). The interior regions of 
the disc (e.g. around i?4 : i) more closely approximate simple 
relative motion between a slowly apsidally precessing disc 
orientation and a rapidly moving orbital frame. 
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Figure 10. Pictorial representation of the mass distribution in the disc for run 9 at four superhump phases. The disc is split into 
100 equal azimuthal sections (dashed lines). Each colour graduation represents 'contours' of particle numbers. Moving outwards, each 
contour represents an increase of 24 particles in each azimuthal section (90 contours in all). Further details are described in the text. 
The outermost dark red points mark the positions of the outermost particles in each azimuthal section. Also plotted are three circles at 
the 3:1, 4:1 and 5:1 resonance radii moving inward respectively (black, dark and light grey) which also act as a guide to the eye to show 
the non-circularity of the contours. 



Figure [12] shows the fitted ellipse parameters for the 
four phases shown in Figure [TO] The average radius of all 
points making up the contour was used. The general trends 
in Figs. [12] and [10] are seen in analogous plots for q — 0.1765, 
including the peaks and troughs in the radial distribution 
of the eccentricity parameter. Within R3.1, the orientation 
of the semi-major axis, ao, changes systematically over the 
superhump period and in the opposite sense to the motion 



of the gas in the disc i.e. precession of the slowly-moving disc 
as viewed from the rapidly rotating binary reference frame. 

Figure [T^] is complementary to Figure [5] Figure \§\ shows 
the eccentricity distributions deduced from instantaneous 
velocities, while Figure [T2l shows the eccentricity distribution 
deduced from the instantaneous mass distribution. Since the 
disc is continuously flexing in a complex way, these are not 
the same. The disc motions can be described as superpo- 
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Figure 12. Fitted ellipse parameters for each 'contour' in Fig- 
ure llOl as a function of radius. See Figure ITul and text for details. 
Results for four superhump phases are presented corresponding 
to the four panels in Figure [10] Vertical lines indicate the 3:1, 4:1 
and 5:1 resonance radii from right to left respectively. Red points 
indicate contours which are not complete. 



sitions of the S(k,l) modes, and the resonance radii act as 
nodes and antinodes in the complex standing wave dynamics 
the disc executes over a full precession period. To summarise 
contributions to these disc motions we compare, in table [TJ 
the strengths of the (1,0) and (2,2) modes. 

For comparison, we show the equivalent of Figures 1101 
and !12l for q = 0.3333 in Figures [TTI and frSl a disc which does 
not show superhumps. Here, as we would expect, the mass 
distribution remains approximately constant over time, and 
the eccentricity is much lower. 



4 PERIOD EXCESS VERSUS MASS RATIO: 
DRAWING TOGETHER OBSERVATION, 
THEORY AND SIMULATION 

If a reliable relationship between t and q can be deduced, 
this would be immensely useful, e can be easily measured 
using relatively modest equipment, while the mass ratio q 
is more fundamental and less easily determined. Patterson 
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Figure 13. As Figure 1121 and corresponding to the panels in 
Figure [TT1 



et al. (2005) fitted observations of eclipsing systems with 
e = 0.18 g + 0.29 g 2 . Figure [141 collates observation, theory 
and simulations of positive superhumps. All observed sys- 
tems with q determined by some means independent of e 
are plotted. Errors in q are formal errors given by the au- 
thors and do not necessarily reflect the uncertainty in the 
method. The eclipsing systems (red circles) should therefore 
be given more weight. We note, however, the scatter of the 
eclipsing systems seems typical of the scatter of the other 
points. Error bars for e denote the range of values observed 
rather than errors in individual values, except in cases where 
only one measurement has been made. Observational data 
is tabulated in Table [S] 

Our high-resolution 3D simulations, which are repre- 
sented by blue squares in Figure 1141 provide a far better 
match with observed systems than previous studies by Mur- 
ray (1998, 2000). As in Murray (2000) we see that simple 
dynamical precession as given by 

3 



■ p(r) 



a, 



(6) 



poorly represents observed systems. This is true even if the 
location of the resonance for a gaseous disc, rather than 
for isolated particles, is used; the location of the resonance 
changes only by < 1 per cent. In a real gaseous disc the 
retrograde effect that pressure forces have on disc preces- 
sion rates must be taken into account (Lubow 1992; Murray 
2000). In a gaseous disc, the excited eccentricity propagates 
through the disc as a wave and is wrapped into a spiral by 
the differential precession of the gas. Murray (2000) assumed 
that the hydrodynamical precession is given by 



+ W v 



(7) 



where ui pr is the pressure contribution to the precession, and 
showed that, under the assumption that the eccentricity is 
tightly wound (that is if it is wound up on a length-scale 
much smaller that the disc radius) then this pressure con- 
tribution at the 3:1 resonance radius is 



2 



Q or b a tan i p 



(8) 
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Table 5. Observed super hump systems with independently determined mass ratio. All periods are given in days. Columns 5 and 6 list 
the minimum and maximum observed superhump periods respectively. Column 7 gives the superhump excess where the errors indicate 
the range of values as calculated from min P sn and max P sn . In column 8, the object type, (E) denotes a WD eclipsing system, while (e) 
denotes a system which shows eclipse of the accretion stream /disc impact region only. 



System 


9 


Pavb 




min P s h 


max P 3h 


e 


Type 


Ref 


OY Car 


0.102(3) 


0.0631209180(2) 


0.06454(2) 


0.064245 


0.06466 


0225+ 0019 
u.uz^a_ 0047 


SU UMa(E) 


1,2,3,4,5 


XZ Eri 


0.1098(17) 


0.061159491(5) 


0.062808(17) 


0.062603 


0.06283 


0270+ 0003 
u.Luru_ 0034 


SU UMa(E) 


6,6,7,7,8 


IY UMa 


0.125(8) 


0.07390897(5) 


0.07588(1) 


0.07558 


0.07599 


0267+ 0015 
u.u-sor_ 0041 


SU UMa(E) 


9, 9, 10 , 11 , 11 


Z Cha 


0.1495(35) 


0.074499 


0.07740 


0.0768 




0.0389_ ooso 


SU UMa(E) 


12 , 13 , 14 , 14 


HT Cas 


0.15(3) 


0.07364720309(7) 


0.076077 






0.0330 


SU UMa(E) 


15 , 16 , 17 


DV UMa 


0.1506(9) 


0.0858526521(14) 


0.08870(8) 


0.0886 


0.08906 


0332+ 0042 

U.VJ3J_ q 012 


SU UMa(E) 


6,6,18,18,18 


OU Vir 


0.175(25) 


0.072706113(5) 


0.07508(9) 


0.07505 




0-0327± 0005 


SU UMa(E) 


19, 19,8,20 


V2051 Oph 


0.19(3) 


0.0624278634(3) 


0.06418(16) 


0.06365 


0.06439 


n 0981 +0-0033 
U - 0JSX - 0.0085 


SU UMa(E) 


21,22,23,23,23 


WZ Sge 


0.060(7) 


0.0566878460(3) 


0.05726(1) 


0.05716 


0.05738 


0101 + 0021 
u.uiui_ 0018 


WZ Sge(e) 


24,25,26,5,27 


VY Aqr 


0.11(2) 


0.06309(4) 


0.06437(9) 


0.0642 


0.06489 


0203+ ' 0082 
u.o/O3_ 0027 


SU UMa 


28,29,5,5,5 


CU Vel 


0.115(5) 


0.0785(2) 


0.08085(3) 


0.0799 




0.0299±j|. 0121 


SU UMa 


30,30,20,31 


SW UMa 


0.14(4) 


0.056815(1) 


0.058182(7) 


0.05790 


0.05833 


0241+ 0026 
0.0050 


SU UMa 


32,33,34,35,36 


HS 2219+1824 


0.19(1) 


0.0599 


0.06184 






0.0324 


SU UMa 


37,37,37 


EK TrA 


0.20(3) 


0.06288(5) 


0.06492(10) 


0.0648 




0-0324± 0019 


SU UMa 


38,38,38,39 


EI Psc 


0.21(2) 


0.044572(2) 


0.04654 


0.04579 




°-0442+2„ la9 


SU UMa 


40,41,42,42 


VW Hyi 


„ ,.,+0.03 


0.074271038(14) 


0.07714(5) 


0.07621 


0.07824 


0.0386± a °ig 


SU UMa 


43,44,44,45,45 


YZ Cnc 


0.22 


0.0868(2) 


0.09204 


0.0905 




°- 0604 t[!.0178 


SU UMa 


46,46,47,47 


WX Hyi 


23+ 07 


0.0748134(2) 


0.07737 




0.0783 


0.0342l;2 0124 


SU UMa 


43,48,49,50 


T Leo 


0.71(15) 


0.0588190(5) 


0.06022(2) 


0.06021 


0.06025 


^ „_„„-+-□ 0005 

0.02381°;™°,! 


SU UMa IP? 


32,51,52 53 52 


U Gem 


0.357(7) 


0.1769061898(30) 


0.20 


0.197 


0.203 


131+ ' 017 
U - 1J1 _0 .017 


U Gem(e) 


54,55,56,56,56 


V603 Aql 


0.22(3) 


0.13809(12) 


0.14640(6) 


0.144854 


0.14686 


0.0602+°;™33 


Fast nova 


57,58, ,59, 60 


UU Aqr 


0.30(7) 


0.163580429(5) 


0.17510(18) 






0.0704 


NL(E) 


61,61,8 


DW UMa 


0.39(12) 


0.136606527(3) 


0.1454(1) 




0.1461 


0.0644;; , 0051 


NL(E) 


62,63,64,63 


MV Lyr 


43+ 19 


0.132335 


0.1377(4) 




0.1487 


0.0405 ±° ° 832 


NL 


65,66,67,68 


AM CVn 


0.18(1) 


0.011906623(3) 


0.012167 


0.012164 


0.012169 


n n9i8+o oo 02 
u.02i»_ 0002 


AM CVn 


69 , 70 , 70 , 70 , 70 


KV UMa 


0.037(7) 


0.1699339(2) 


0.170529(6) 


0.17049 


0.17073 


0035+ 0012 

u,UUM -0 .0002 


BHXRT 


71,72,73,74,8 


(XTE J1118+480) 

QZ Vul 0.042(12) 


0.3440915(9) 


0.3469(1) 




0.3474 


0.0082l;2 0014 


BHXRT 


75 , 75 , 76 , 77 


(GS 2000+2) 
V1487 Aqr 


0.058(33) 


30.8(2) 


31.4 


31.2 


31.6 


0195+ ' 0065 
u.ui»o_ 0065 


BHXRT 


78 , 79 , 79 , 79 , 79 


(GRS 1915+105) 


















V518 Per 

(GRO J0422+32) 

GU Mus 


„ 111 + 0.027 
" ,iu -0 .033 

0.13(2) 


0.2121600(2) 
0.432602(1) 


0.2157(10) 
0.4376(10) 






0.0167 
0.0116 


BHXRT 
BHXRT 


80,80,81 
82,83,77 


(N Mus 1991) 
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Figure 14. Superhump period excess plotted as a function of binary mass ratio for both observed systems and for SPH simulation. Also 
plotted are the dynamical and hydrodynamical theoretical predictions. Simulations presented in this paper are displayed as filled squares, 
as detailed in the top right legend. The second legend down refers to previous works by Murray. The third legend refers to theoretical 
predictions, and the legend for observational data is in the bottom right. The inset shows data over a large range of q, whilst the main 
panel focuses on the data at low q where most of the points are clustered, and is plotted on a logarithmic rc-axis. 



where i p is the pitch angle of the spiral wave, and c s is the 
sound speed. 

For each eclipsing SU UMa system we have calculated 
an inferred pressure contribution to the precession rate, ui pT 
(column 7 in Table [3]l in the manner of Murray (2000). A 
weighted mean gives lo pi — — 1.36 radd" 1 . This is combined 
with the mean mass of the WD in systems below the period 
gap (Smith & Dhillon 1998) to calculate predicted hydro- 
dynamical precession rates assuming that the precessional 
pressure contribution is similar for all systems. This is plot- 
ted as a solid curve on Figure 1141 We note that Pearson 
(2006) considers the inclusion of pressure effects in an al- 
ternative way. Whilst a reasonable fit to the observations is 
achieved, the fact that there is a distribution in the eclipsing 
systems above and below this curve, demonstrates that the 
situation is not so simple. Possible contributory factors are 
distributions in primary mass or in disc temperature. We 
have plotted two further curves on Figure PHI one represent- 
ing the hydrodynamical prediction for larger WD masses, 
and the other encompassing both larger Mi and higher ac- 
cretion disc temperature. These are the mean value of Mi 
for the eight SU UMa eclipsing systems and the value of ui pl 



found by Murray (2000) respectively. We tabulate Mi, the 
total system mass, Mt, and the mid-plane sound speed at 
the resonance radius for each of the eclipsing SUUMa sys- 
tems (Table [3}. These numbers do not seem to provide any 
clue as to the true cause of the spread in observed systems. 
In calculating c s , we assumed the same mass transfer rate 
throughout (10 -9 Mq yr _1 ). Undoubtedly this, and conse- 
quently the temperature, varies from system to system. 

Figure [T3] also shows the best fit that Goodchild & 
Ogilvie (2006) found to their analytic curve. They found 
that retrograde pressure terms are in fact only ~ 1 per cent 
of the dynamical term, far smaller than the values we infer 
above. They suggest that the offset of observation from dy- 
namical precession rate be due instead to averaging over the 
disc, the eccentricity being distributed throughout the disc 
rather than being sharply peaked at the resonance itself. 
This is in accordance with our findings that dissipation- 
powered superhump overwhelmingly originates in the disc 
regions within R3.1. Goodchild & Ogilvie (2006) found that 
their eccentricity distributions peaked close to 0.37 a, a value 
put forward by Patterson (2001) to match observation. We 
have plotted the dynamical curve as evaluated at 0.37 a (dot- 
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dashed light-blue line). We have also plotted two further 
lines which show dynamical curves evaluated at some frac- 
tion of the resonance radius, one to give the best weighted 
fit to the eclipsing SUUMa systems (red) and the second 
to give the best fit to our 3D simulations that have reached 
equilibrium, runs 5 to 10, (blue). We find values of 0.83J?3 ; i 
and 0.87i?3:i respectively. 



5 DISCUSSION 

The precession rates for the simulated discs provide a much 
closer match with observation than has been achieved previ- 
ously. There are a number of reasons for this. The improved 
mass resolution changes the stream-disc impact, leading to 
a different angular momentum distribution in the disc and a 
different precession rate. The radial sound speed distribution 
is correct for a steady-state disc (as opposed to an isother- 
mal one as in the case of Murray (1998)), which means 
that the variation of density with radius is more realistic, 
which in turn will determine whether the disc precesses and 
at what rate. These updated simulations have hotter discs. 
This means the retrograde effect of pressure on precession 
will be greater (Lubow 1991a). The viscosity is more in line 
with what is inferred for an a-disc in the high state. The 
simulations of Murray (1998) were too viscous, which al- 
lowed the disc to grow and penetrate the resonance too eas- 
ily. The propagation of the eccentricity inward through the 
disc would also have been inhibited. If the precession rate 
is dictated by a weighted average of the eccentricity as sug- 
gested by Goodchild & Ogilvie (2006) , then the inhibition of 
inward eccentricity propagation would lead to high preces- 
sion rates. Finally, the extension to 3D changes the character 
of the resonance. We can see in Figure |S] that the accretion 
disc look very different in 2D and 3D. One reason for this is 
that the character of the stream-disc interaction differs: for 
2D mass injection all the particles follow one another exactly 
and so the stream tends to punch through the outer disc and 
deposit its angular momentum at r < R3.1. Conversely in 
the 3D case the particles are more easily captured by the 
outer disc and so they effectively reduce the specific angu- 
lar momentum of the outer disc. The growth and decline of 
the eccentric mode is strongly influenced by the interaction 
with the mass transfer stream. The picture is further com- 
plicated by the results of Kunze, Speith & Hessman (2001) 
who found substantial stream overflow in their simulations. 
An important extension of this work would be to system- 
atically isolate the importance of each of the effects which 
contribute to the disc precession rate. 

It appears that the offset of both observation and simu- 
lation from dynamical expectations of precession rates is due 
to the averaging over the radii interior to R3.1 which partici- 
pate in the disc precession, and the inclusion of a retrograde 
effect of pressure forces is insufficient (Goodchild & Ogilvie 
2006). In Section [3] and Figure [5] we show that the eccen- 
tric instability is indeed manifest at radii as small as 0.15a, 
so the majority of the disc area contributes to the mean 
precession rate. We note that Goodchild & Ogilvie (2006) 
required a very low semi-thickness parameter (h = 0.003) to 
fit observations, whereas in our simulations this is an order 
of magnitude higher (Table [2)| and in better agreement with 



observational constraints on, and theoretical expectations 
of, disc thickness. 

Empirically, systems which show superhumps generally 
have g < 0.24. There are two classes of exceptions: nova-likes 
and U Gem, and magnetic systems. For the former, Osaki 
(2005) pointed out that Paczynski's derivation of iitidea as- 
sumes the disc to be cold. A sufficiently hot accretion disc 
may expand beyond this owing to a weakening of shocks 
by strong pressure effects at the last non-intersecting orbit. 
Osaki (2005) argues the persistently high-state nova-likes, 
and the unusually long 1985 outburst of U Gem in which su- 
perhumps were reported, may satisfy this temperature cri- 
terion. The magnetic systems, TV Col and possibly TLeo, 
may have discs which are pushed out by magnetic forces 
(Retter et al. 2003). 

The high resolution 3D simulations presented here re- 
produce this distribution well; the upper limit for which su- 
perhumps are observed is also q ~ 0.24, albeit at this mass 
ratio it is a very slow process indeed. This also compares well 
with the theoretical upper limit of q < 0.25. For q — 0.2422 it 
appears that the eccentricity starts to grow more than once, 
but fails each time and the eccentricity falls away again, 
presumably because the encounter with the resonance was 
only marginal. The final eccentricity is highest for q ~ 0.1, 
and declines gradually as q becomes less extreme. This is 
expected as the disc has more room to grow beyond the 3:1 
resonance before being tidally truncated for extreme mass 
ratios. For the two most extreme mass ratio simulations, 
with g = 0.0526 and g = 0.0256 the final eccentricity is neg- 
ligible. In the q — 0.0256 case there is no eccentricity growth. 
The q — 0.0526 case is interesting: as the disc is growing, the 
eccentricity initially increases and superhumps are appar- 
ent, albeit weak and ill-formed (Figure |3). This eccentricity 
is, however, damped away again. What is the source of the 
damping? One possibility could be the action of the 2:1 reso- 
nance which may be excited in ultra-low mass ratio systems 
and which can act to damp eccentricity (Lubow 1991b). Per- 
haps this results in a close competition between the 2:1 and 
3:1 resonance where the 2:1 resonance is marginally excited 
in the case of q = 0.0526, and the 2:1 resonance only comes 
out on top once the disc has grown and this resonance is 
sufficiently populated. In reality, the most extreme mass ra- 
tio system in which superhumps are observed is the black 
hole X-ray transient KVUMa (XTE J1118 +480), which has 
q = 0.037 ± 0.007. In 2D, the range for which the simulated 
discs become eccentric extends to much higher mass ratios 
(g «5 0.4815) (Table □ Figure [HJ). The confinement to 2D 
means that the character of the resonance is different and 
the strength of the resonance is increased. 

We find that the growth rate of the eccentricity (the 
growth rate of Snm) i s highly dependent on g, with high 
mass ratio systems taking a very long time indeed to be- 
come eccentric (Table [TJ. We see in Figure [2] that our new 
simulations have slower growth rates than previous 2D sim- 
ulations by Murray (1998). Observationally, the rise times 
of superoutbursts is of order a day or so, and superhumps 
are generally detected within a day or so. The outburst on- 
set is, however, probably governed by the thermal-viscous 
instability and hence its timescale is independent of the disc 
eccentricity. Generally there are no observations suitable for 
assessing whether or not the disc was eccentric before the 
superoutburst began, with intensive observations beginning 
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after the rise to outburst, so timescales on which SU UMa 
discs typically develop their eccentricity is ill-constrained. 
We note further that our simulations do not include the 
thermal-viscous instability, so their evolution differs from 
that of SUUMa discs. Lubow (1991a) found an analytical 
expression for the growth rate of the eccentricity, finding it 
to be proportional to q 2 . This is applied to an ideal narrow 
fluid ring. Clearly the eccentricity growth rate of our simu- 
lated discs is not proportional to q 2 . Osaki (2005) uses the 
dependence of eccentricity growth on q 2 to propose, as a re- 
finement to the TTI model, an explanation for type A/type 
B superoutbursts, namely those which show a precursor and 
those which do not. He suggests that it depends on whether 
the eccentricity growth rate is large enough to excite a signif- 
icantly eccentric disc within the duration of a normal out- 
burst, and that this is why most SUUMa systems having 
relatively low mass ratio show only Type B superoutbursts. 
However, if it is the case that high q means lower growth 
rates as we see, then this argument fails, and probably other 
factors contribute. 

Goodchild & Ogilvie (2006) also find extremely low 
growth rates in their analytical work. They use the growth 
term found by Lubow to describe the rate at which eccentric- 
ity is created at the resonance, but consider further how this 
eccentricity propagates through the disc. They explain their 
low growth rates as due to the eccentricity being strongly 
suppressed at the resonance itself. They find that the growth 
rates depend most strongly on mass ratio and on bulk viscos- 
ity, with further weaker dependence on disc semi-thickness. 
In column 7 of Table [5] we list the parameter, h, referred 
to by Goodchild & Ogilvie (2006) as the characteristic disc 
semi-thickness, and in column 11 the bulk viscosity at R3.1 in 
our simulations. For the 3D simulations, these are ~ 0.036 
and 0.2 respectively. Comparing our Figure [2] with Good- 
child & Ogilvie (2006) 's results, we see that the trend in our 
3D points matches quite well with the h — 0.01 line (the 
highest value of h given) in their figure 10. In particular we 
see in both cases low growth rates at high q, a maximum 
at q ~ 0.08 and a steep drop in growth rate at mass ratios 
below this. Our growth rates, though, are about a factor of 
10 higher. Looking to their figure 9, this could be explained 
by our high semi-thickness parameter. Our bulk viscosity, 
though, is also rather high. It would be very interesting to 
make a study of empirical growth rates and their depen- 
dence on q and on other known parameters to compare with 
these findings. This would require systematic monitoring of 
dwarf novae to catch the onset of outburst. 

There are, however, distinct differences between the 
work of Goodchild & Ogilvie (2006) and the simulations 
that we present. Tidal modes, which would presumably act 
to truncate the eccentric mode, are not included in their 
work. The outer boundary conditions differ. We see too in 
Section [3] that the eccentricity distribution as a function of 
radius appears quite different from their findings. The situa- 
tion in the simulations is further complicated by the presence 
and importance of the tidal 2-armed spiral structure which 
is not included in Goodchild & Ogilvie (2006) 's work. 

In our simulations we see the superhump period de- 
creasing (Figure [TJ as is often observed over the course of 
a superout burst (Patterson et al. 1993a). As the period 
changes, the disc eccentricity is increasing, and, in most 
cases, the disc mass has begun to decrease in response to 



the enhanced tidal torques. The superhump period decrease 
in the simulations can then be explained by the eccentric 
wave propagating inward, and additionally by radial shrink- 
ing of the disc. We are unsure how to explain the exception, 
q — 0.0526, which shows an opposite behaviour for the ini- 
tial part of the simulation. Perhaps it is related to the ideas 
of Uemura et al. (2005) , where they suggest an explanation 
for +ve Pah observed in a few cases. They suggest P s h is 
related to the amount of matter beyond R3.1, so allowing 
for an outward propagation of the eccentric wave. For low 
q, the distance between R3.1 and -Rtidcs is greater. 

The observational data points to a many-valued t(q) 
relation (Figure I14p . In particular 3 of the BHXRTs show 
systematically lower precession rates than those of CVs. 
QZ Vul (GS 2000+2) might be an exception to this. However 
the superhump period measurement is uncertain as it has 
been sparsely observed. The microquasar V1487 Aqr (GRS 
1915+105) is a much longer period system and may not 
be comparable to the other BHXRTs. The accretion discs 
in BHXRTs are irradiated by the central X-ray source. We 
would expect these discs to be both hotter and also thicker 
due to the bloating effect of irradiation. Both of these would 
act to reduce the precession rate, due to the retrograde ef- 
fect of pressure and due to the dependence on semi-thickness 
found by Goodchild & Ogilvie (2006). It would be very inter- 
esting if e and q for further BHXRTs could be determined. 
The only ultra-compact helium binary included, AMCVn, 
also lies below the main cluster of points in Figure 1141 As 
Roelofs et al. (2006) noted, the helium accretion disc in 
AM CVn could be thicker than its hydrogen-rich counter- 
parts. 



6 SUMMARY 

The main findings in this work can be summarised as follows: 

• We present improved accretion disc simulations for a 
range of q. The main improvements are both numerical and 
physical: a higher mass resolution, extension to 3D, more 
realistic disc temperature and viscosity, and a radial depen- 
dence of sound speed appropriate to a steady-state accretion 
disc. We ran the simulations until equilibrium was reached. 
For 0.08 < q < 0.24 the 3D discs reach an eccentric equilib- 
rium and show a superhump signal in their energy dissipa- 
tion rate (which we refer to as a simulated lightcurve). 

• The e(q) dependence for the SPH simulations presented 
in this work shows a greatly improved match with observa- 
tion than previous simulations. 

• No high resolution 3D disc with q > 0.24 developed 
superhumps. This agrees with theoretical expectations and 
matches the majority of observations. 

• The region of the disc within R3.1 is overwhelmingly re- 
sponsible for generating the dissipation-powered superhump 
in the simulations. 

• If the difference between observed precession rates and 
dynamical precession rates calculated at the 3:1 resonance 
radius is due to averaging over the disc as Goodchild & 
Ogilvie (2006) suggest, then we find that the best fit char- 
acteristic radius of the eccentricity distribution at which the 
dynamical precession rate is evaluated to be 0.87-R3 : i and 
0.83i?3:i for the 3D simulated discs and the observed eclips- 
ing systems respectively. The differences between these two 
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best-fit radii may be partly due to the differing surface den- 
sity distributions in the two cases. 

• Our simulations show the effect of the increased effi- 
ciency of tidal return of angular momentum to the binary 
for an accretion disc which has become eccentric. The disc 
mass approaches a new lower steady-state value as the disc 
becomes eccentric. This is exactly as asserted by the TTI 
model. With the assumption of a radial dependence of vis- 
cosity, we deduce an effective ~ 4 per cent increase viscous 
torque between a disc which is circular and one that is ec- 
centric. The increase depends on q. 

• As the eccentricity grows and the disc mass falls, the 
superhump period decreases. 

• The dependence of eccentricity growth rates on q that 
we see in the simulations presented here is comparable to 
the work of Goodchild & Ogilvie (2006). Particularly, we 
find that for high mass ratios the growth rates are very low 
indeed, in contrast to the result of Lubow (1991a). This 
needs to be reconciled with observation. 

• We show that superhumping discs have noticeable ec- 
centricity even in their inner regions (r ~ 0.15a). Conversely, 
non-superhumping discs are seen to be eccentric only in 
their outer regions. In this case however, this 'eccentricity' is 
steady-state and has origin in tidal distortions, being there- 
fore different from that which dominates the main body of 
the superhumping discs. We characterise the eccentricity dis- 
tributions using two different methods. 

• The disc motions can be described as superpositions of 
the S(k,l) modes, and the resonance radii act as nodes and 
antinodes in the complex standing wave dynamics the disc 
executes over a full precession period. We characterise the 
disc motions on P a h, the key timescale for the powering of 
the observed superhumps. 

• The 4:1 and 5:1 resonances may play roles in the dy- 
namics of eccentric discs for g < 0.24. This may explain 
why the observed precession rates are closer to the dynamic 
precession rate at the 4:1 resonance than they are to the 
dynamic precession rate at the 3:1 resonance. 

• The observational data shows a multi-valued e(g) rela- 
tion. In particular, the BHXRTs show systematically lower 
precession rates than those of the CVs, which may be ex- 
pected when the higher temperature and thickness of their 
irradiated discs is considered. 
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